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ABSTRACT 


The work presented here is a further development of the ideas ex- 
pressed in this author's earlier work on the modification of Hansen's 
lunar theory. We introduce here several devices which permit us to 
reduce the number of quantities which must be re-computed with each 
iteration cycle. We refer the motion of the satellite to an "almost ideal" 
system of coordinates with two axes lying in the mean orbital plane. 
These axes are selected in such a way that the differential equations of 
the relative motion differ from the equation of motion in Hansen's ideal 
system only by the terms of higher orders. The position of the satellite 
relative to the almost ideal system is determined by Hansen's coordi- 
nates l + v and n 0 Sz of the projection of the satellite on its mean or- 
bital plane and by the "elevation" £ of the satellite relative to this 
plane. The position of the almost ideal system itself with respect to an 
inertial system is determined by two uniform rotations, one around the 
normal to the orbital plane of the sun, and one around the normal to the 
mean orbital plane of the satellite. The problem of integration is re- 
duced to solving a linear partial differential equation by means of suc- 
cessive approximation, or to expanding the integrating operator into a 
series of products of two linear operators. One of these operators is a 
linear partial differential operator and the other is the inverse of a 
linear partial differential operator with the constant coefficients. We 
propose here, as in the previous work, to compute Hansen elements 3, 
T and V separately and to fuse them into w bypassing Hansen's 
function W. 

The application of Hansen's theory to the Jovian satellites is a topic 
of considerable astronomical interest and importance. M. Charnow has 
programmed the modified Hansen's theory in order to produce the 
yearly epheme rides of the X th satellite of Jupiter. 
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ON SOME POSSIBLE SIMPLIFICATIONS AND 
CHANGES IN HANSEN’S LUNAR THEORY 


by 

Peter Musen 

Goddard Space Flight Center 


INTRODUCTION 

This work represents a further development of possible modifications of Hansen's (1838, 1862) 
lunar theory which were presented by the author in a previous work (Musen, 1963). We direct our 
attention to Hansen's lunar theory because it is applicable with some changes to satellites moving 
in orbits with high inclinations relative to the orbital plane of the Sun. In particular, its possible 
application to the motions of Jovian satellites (Mullholland, 1965; Charnow, 1967) is a topic of con- 
siderable astronomical interest and importance. 

The main kinematical characteristics of Hansen's theory are: 1) the separation of the periodic 
perturbations in the orbital plane from those of the orbital plane itself; and 2) the use of a fixed 
rotating ellipse as an intermediate orbit. Normally, the periodic perturbations in the orbital plane 
are larger than the periodic perturbations of the orbital plane itself. These two types of effects are 
not completely independent, but their mutual influence is of higher order and it can be calculated. 

The differential equations of motions of the satellite in its osculating orbital plane are referred 
to a rotating reference frame rigidly connected to this plane. This frame of reference Hansen terms 
as "ideal," because the relative and the absolute velocities coincide and because the Coriolis and 
the centrifugal forces cancel each other. Thus, the vectorial differential equation of the relative 
motion has the same form as in an inertial system. An intermediate orbit is an ellipse of constant 
shape located in the osculating orbital plane and rotating with constant angular velocity relative to 
the ideal frame around the axis normal to this plane. A fictitious "auxiliary" satellite moves in 
this ellipse in accordance with Kepler's laws with respect to pseudo-time. The position vector of 
the real satellite at any given time is determined as its deviation from the position vector of the 
auxiliary satellite in time and space. 

In Hansen's theory all the angular perturbations in the orbital plane are combined into one 
single angle n 0 8z, the perturbations of the mean anomaly. The radial perturbations are given by 
1 + v, the ratio of the radius vector of the real satellite to the radius vector of the fictitious satel- 
lite. The position of the osculating orbital plane and, consequently, of the ideal system of coordi- 
nates is determined by the mean position of the orbital plane and by three additional parameters 
which absorb all the small periodic oscillations of the osculating orbital plane around its mean 
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position. The final output in the classical form of Hansen's theory are the trigonometric expan- 
sions of the perturbations in the sine of the latitudej of 1 + v and of n 0 Sz. 

The author has suggested (Musen, 1959) that the expansion of the sine of the latitude be dis- 
carded if the inclination is large. Instead, he proposes to use a set of Euler’s four "redundant” 
parameters to determine the periodic deviations of the osculating orbital plane from its mean po- 
sition. Euler's parameters, like Hansen's, carry all the periodic effects in the node, the inclina- 
tion, and the position of the departure point. The position of the mean orbital plane itself relative 
to the inertial frame is influenced only by the secular effects in the node, the perigee and the de- 
parture point. The use of Euler's parameters makes all the angular arguments of the theory linear 
relative to time from the outset and the expansion of the disturbing function becomes more sym- 
metrical, as well as more algebraic in form, as compared to the expansion in terms of Hansen's 
parameters. 

Developing the theory of satellites in terms of Brendel's coordinates (Brendel, 1925), the 
author has noticed (Musen, 1967) the existence of an "almost ideal" reference frame in the mean 
orbital plane, such that the differential equations of motion of the satellite relative to this frame 
differ from the equations in the ideal system only by the terms of higher orders. The use of the 
mean orbital plane instead of the osculating one, as well as the introduction of the "almost ideal" 
system, permits one to discard the use of Hansen's or Euler's parameters. In addition no "re- 
dundant" parameters appear in the reformulation of the theory. The position of the satellite is 
determined by Hansen's coordinates of its projection on the mean orbital plane and by its "elevation” 

£ relative to this plane. This elevation is very small and its use permits one to contract the iteration 
process considerably. This system leads to an expansion of the disturbing function such that the 
coefficients of £ and 1 + v remain unaltered from the iteration to iteration. This circumstance, 
consequently, represents a substantial simplification as compared to the system based on the use 
of Hansen or Euler parameters. 

Hansen's original lunar theory contains the derivatives of the disturbing function relative to 
the eccentricity and relative to the parameters associated with the position of the osculating orbital 
plane. Evidently, the literal expansions of the disturbing function and of its derivatives are required 
before the process can be made completely numerical. We retain here our previous system (Musen, 
1963) which makes use of expansions of the components of disturbing force and which somewhat 
resembles Hansen's planetary theory. This system permits the numerical substitution of the mean 
eccentricity and of the mean inclination from the very start. We expand the disturbing force into 
an infinite series in terms of the mean anomalies of both the satellite and of the sun, and in powers 
of 1 + v and £ . Expansion in terms of the disturbed eccentric anomaly E of the satellite and of the 
true anomaly f ' of the sun is also possible. The advantage of this last expansion is that the dis- 
turbing function can be represented as a series of trigonometric polynomials in E and f ' . However, 
if the higher order effects are being sought, then the problem of integration is much simpler if the 
mean anomalies are employed. The determination of the mean motion of the perigee in Hansen's 
theory involves division by the eccentricity. Thus, an inconvenience arises if the mean eccentricity 
is very small. The situation can be remedied by keeping the first power of the eccentricity in 
literal form— in other words, by representing each term of the trigonometric series as a "vector” 


2 



Ae + B, with the following laws of multiplication for the "unit vectors": 


e ■ 1 - 1 • e - e , e 2 - ee Q , 

where e 0 is the numerical value of the mean eccentricity. The motion of the perigee is obtained 
from the condition that no constant for the form eA is present in the derivative of Hansen's element 
>P. After the work is completed e should be replaced by e 0 . 

The main charm of Hansen's lunar theory is that all the perturbations in the orbital plane are 
determined by one single function W. As an auxiliary step Hansen introduces a second function w. 

In this second function the "elliptic" mean anomaly, which describes the motion of the fictitious 
satellite, and the mean anomaly taken as the argument in the expansion of perturbations are sep- 
arated. Then the differential equation for w depends only upon the derivatives of the osculating 
elements, but not upon the motion of the fictitious satellite directly. After the integration of the 
expansion of dw/dt is completed, the distinction between the two types of mean anomalies is removed; 
and by applying Hansen's "bar-operator" w becomes w. With respect to these two functions we take 
the same approach as before (Musen, 1963): only W has a direct kinematical meaning ; w is merely 
an artificial device to blend three independent series into one. 

However, neither the number of terms nor the computing time can be diminished by using w 
instead of the independent three series. We propose here, as in the previous work, to compute the 
three series for Hansen's elements 5, T, t separately, without forming the W-function. After they are 
computed, we fuse them together to form the w-function. We suggest the method of iteration to 
solve the problem, partly because the starting approximation might include some perturbation ef- 
fects already, and thus the input information might go beyond Hansen's moving ellipse, and partly 
because the programming is more uniform if the iteration process is employed. An additional 
reason for employing the process of iteration is its close association with the fixed point theorems 
in the modern theory of differential equations. These theorems are now being used widely to es- 
tablish the existence and the properties of solutions. 

A considerable portion of Hansen's lunar theory constitutes the expansion in powers of the 
perturbations in the mean anomaly n 0 Sz . Almost every step requires this expansion and the de- 
termination of w even requires the expansion in powers of n 0 Sz twice. We can abbreviate the proc- 
ess of expansion considerably by developing the theory in terms of the disturbed mean anomaly 
directly without resorting to the expansion in powers of n 0 Sz at almost every step. Then the problem 
of integration relative to time can be reduced either to solving of a linear partial differential equa- 
tion by means of successive approximations or by expanding the integrating operator (d/dt )“ 1 into 
a series of products of two linear operators. One of these operators is a linear partial differential 
operator and the other is the inverse of a linear partial differential operator with constant coef- 
ficients at the operators of the partial derivatives with respect to the four basic arguments. 

Several numerical inaccuracies were found in Hansen's expansions of the coordinates of the 
Moon. However, this circumstance does not affect the correctness of his basic theoretical thought. 
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The theory is applicable to the natural planetary satellites when solar effects are dominant. Several 
small, but observable and physically important effects of non-solar origin, which are so important 
in the motion of the Moon, have no significance in the computation of yearly ephemerides of the 
satellites of the outer planets. 

The Hill-Brown lunar theory makes use of the variational solution as an intermediary orbit. This 
solution is free from the influence of the eccentricity and of the inclination and represents a natural 
start if one wants to expand the coordinates and the frequencies in powers of the constants of in- 
tegration. This does not mean, however, that in general the variational solution contains the effects 
which are the most important numerically. The choice of the rotating ellipse as an intermediary is 
a natural one in a purely numerical approach, because this ellipse introduces the mean motions of 
the node and of the perigee from the start, and because the large terms in the perturbations have 
the eccentricity of the satellite as a factor. 

The unorthodox manner of treating the perturbations, the absence of the direct appeal to the 
method of the variation of the astronomical constants and, especially, the use of the "disturbed" 
time, and the rather difficult style of Hansen's exposition, make the reading of Hansen's work rather 
demanding to modern readers. But the reader willing to overcome these obstacles will in the final 
analysis find Hansen's basic idea to be admirably clear and well adaptable to the use of modern 
electronic machines. M. Charnow (1966) has programmed an accurate expansion of Hansen's per- 
turbations using the modified Hansen's theory as developed by the author. On the basis of his ex- 
pansion the yearly ephemerides of the X th satellite of Jupiter are being computed. Since 1958 this 
satellite was not observed because of the absence of the accurate theory and of the ephemerides. 

On the basis of Charnow's ephemeris the satellite was re-found in 1967 by E. Roemer in Arizona. 

The application of this expansion to the XI th satellite of Jupiter is being planned. 

DIFFERENTIAL EQUATIONS OF MOTION OF THE SATELLITE 

Let k be the unit-vector normal to the orbital plane of the planet and let us consider this plane 
as fixed. Let P be the unit-vector along the ascending node of the mean orbital plane, R be the 
normal to this plane and Q = RxP. The mean elements are selected in such a way that no purely 
secular or mixed terms will appear in the expansions of the rectangular coordinates. The mean 
orbital plane rotates around k with constant angular velocity 


n 0 r - n 0 yk 

where n 0 , the mean anomalistic motion of the projection of the satellite on the mean orbital plane, 
must be given in advance as a constant of the theory. Let I be the inclination of the mean orbital 
plane toward the orbital plane of the planet; it is also one of the basic constants of the theory. We 
have: 


n 0 ? = n 0 /3Q + n 0 aR 


(1) 
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where 


a = y cos I , P ' 7 sin I. (2) 

Let E be the mass of the planet, m' the mass of the sun, r the position vector of the satellite, and r' 
the position vector of the sun. The gravitational constant is put equal to one. We consider the mass 
of the satellite as negligible. The mutual distance | r ' - r | between the sun and the satellite we des- 
ignate by A . The equation of motion of the satellite relative to the rotating frame (P, Q, R) has the 
form: 


• d 2 r dr 

■^2 + 2n orx at 



+ E grad r R 


(3) 


where 


R 


m ' / 1 r • r ' \ 

-E[K- r ,T ) 


(yx r ) ‘ 


(4) 


Let p be the projection of r on the mean orbital plane and £ be the "elevation" of the satellite 
with respect to this plane. We have 


r = P + i R (5) 

Substituting Equation 5 into Equation 3 and taking Equation 1 into consideration, we have 


d ^ + 

dt 2 


dp d£, 

2 n o n R * dt + 2n o at 



E grad p R 


(6) 


and 


d 2 4 
dt 2 


2n 0 fiP 


dp 

dt 



^ 1311 
E ‘5T 


(7) 


In the mean orbital plane we choose a second reference frame rotating with the constant velocity 
-n 0 a relative to the frame (P, Q) around the R -vector. In accordance with the basic rules of kine- 
matics we have to replace in the Equations 6 and 7 


dp 

dt 


by 


dp 

dt 


n Q aR x p , 
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d 2 P 
dt 2 


, d 2 p dp 2 2 

by -^- 2n o aR>< aF~ n o a / > ' 

where dp/ At and d 2 p/ dt 2 now designate the relative velocity and the relative acceleration of the 
satellite in the new rotating system. 

Equation 6 takes the form: 


d 2 p 
dt 2 


E 

— r p + E$ , 

P 


where 


2 n o a „ d£ 

* - grad^n - g— p dt tan I 


( 8 ) 


0) 


and 


n 


m 1 
~E 



a 2 p 2 ] . 


( 10 ) 


Let v, t be the coordinates of the satellite relative to the system (P, Q, R): then 

(ID 
(12) 

Let us define the mean semi-major axis of the orbit of the satellite by means of the equation 


r - dfP + 7)Q+£R, 

p = <f P + tjQ . 


( 13 ) 


We have also for the corresponding equation for the sun 


n 


' 2 
o 


a 


/ 3 

0 


m* + E 


( 14 ) 


or, in a more convenient form, 


l - 


m' + E 


a 


' 3 
0 


6 



Substituting Equations 11 through 14 into Equations 9 and 10 we deduce: 


2a d£ 

a o* ~ & rad p a o n i P dt tanI 

n o a o 


( 15 ) 


and 


a 0 n 




n 

i\ 

a 2 ri 

It 

pi 

a 2 t-2 

a o 


(16) 


where 


and a 0 is the parallactic factor 


m 



All the terms in Equations 15 and 16, excepting the first term in fi, are of higher orders. Thus the 
equation of motion in the form given by Equation 8 differs from the corresponding equation of mo- 
tion in a nearly ideal system by the terms of higher orders. This is why we designate the second 
rotating system of coordinates as an "almost ideal." 

Equation 7 can be written in the form: 


d 2 C, dO. d£ 

~ E Jl + 2n 0 a gj tan I . 

Making use of Equation 13 we can write this equation in the form: 

d2 (£ /a o) „ 5 ( a 0 n ) d ^ /a o 

dt 2 = n 0 djW^] + 2n 0 “ dt tanI ' 


(17) 
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To facilitate the integration we put this last equation in the form: 


d2 U /a o) E i 9 ilfoV) d ^ /a o 

~ It’ = 'T 3 ' ^ + "o ^7^) + 2 n 0 a tan I 


( 18 ) 


where 




a„n + 


1 /M 3 av 

2 \ p / V a o/ 


The Equation 8 of motion relative to the almost ideal system of coordinates of the projection 
of the satellite on its mean orbital plane has the same external form as the equation of the disturbed 
planar motion in an inertial or in an ideal system. It is clear, that because of this identity of forms 
the process of integration of Equation 8 will go formally along the same line as, say, in the ideal 
system of coordinates. Similarly, as in the ideal system, all theory of osculation, and the defini- 
tions of the elements as well as of Hansen's intermediary ellipse and of Hansen's coordinates, can 
be transferred automatically to the motion in an almost ideal system. Of course, the elements 
which we introduce are rather the constants of integration and they are not identical to the osculat- 
ing elements of motion in an inertial system. In the following exposition when referring to the 
"elements" we shall understand consistently such "relative" elements in the almost ideal frame. 

The intermediary orbit in our theory is the ellipse of constant shape rotating uniformly with 
constant angular velocity n 0 y relative to the almost ideal system around the normal to the mean 
orbital plane. The auxiliary satellite P t is moving in this ellipse in accordance with Kepler's laws 
and its position in the ellipse is determined by the standard set of equations: 


— cos f - cos E - e n , 
a o u 

(19) 

~z sin f = Vl - e„ 2 sin E , 
a o f 0 

(20) 

- 1 - e 0 cos E , 

a o 0 

(21) 

E - e Q sinE = n Q z(t) + c 0 , 

(22) 


where a 0 is the semi-major axis of the auxiliary ellipse, e 0 its eccentricity, and r, 7, E, n 0 z + c 0 
respectively the radius-vector, the true, the eccentric and the mean anomalies of P x . At a given 
time t the projection of the real satellite P on the mean orbital plane will have the same direction 
as the position of P t at the time z. Designating by r(t) the position vector of p, by p(t) the position 
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vector of the projection of P on the mean orbital plane, and by r(z) the position vector of P, , we 
can write the basic relations of the theory in the form: 


r(t) 

- p( t ) + ? R , 

(23) 

PC*) 

= (1 + V) p(z) , 

(24) 

Z 

= z(t) - 

(25) 


The position of the perigee of the auxiliary ellipse relative to the x-axis of the almost ideal system 
is 77 0 +n 0 y t , where n 0 is a constant. Designating by v the polar angle of the projection of P with 
respect to the x-axis of the almost ideal system, we can split Equation 24 into two standard rela- 
tions of Hansen's theory: 


p = (1 + v) p (26) 

v r 7 + 77 0 + n 0 yt . (27) 

We introduce now the elements osculating relative to the almost ideal system in a sense we 
discussed before. Let a be the relative osculating semi-major axis, e the relative osculating ec- 
centricity, x the longitude of the relative osculating pericenter reckoned from the x-axis of the 
almost ideal system. The osculating mean motion is defined as usually by the equation: 


E 

a 3/2 


We also make use of Hansen's elements 


0 

h 



(28) 


(29) 


and 


1 h + 2 h 0 3e o h 0 


ecos (x' 7r 0 ' n o5 ,t ) ~ e o 


1 - e 2 


*0 h 3 
1 h + 2 h 0 2 e o T ’ 


e cos (y “ "o ~ n o y 1 ) - e o 


!- e o 2 


(30) 

(31) 
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( 32 ) 


¥ 



e sm * 


yt) 


where 


and, similarly, 


h = — ^ 

^ a ( 1_ e2 ) 


n a 

"^1 - e 2 


rfj 

/ a o ( 1 _ e o 2 ) 



(33) 


(34) 


Designating by ES the radial component of the disturbing force and by ET its orthogonal component 
we have (Hansen, 1838) with some small changes in notation: 


d hp 
dt h 



(35) 


dT 

dt 



1 - e 2 
1 e o 


cos f + 


1 _e o 2 


K 


P K t ) + T sin f ‘ p ( a o s )[ 


(36) 


dt 


h 

T+ 2 TT 

1 1 n 



sin f ■ p[ a„ T) - cos f • p(a Q s) 


(37) 


From the expression (15) for a 0 * we obtain: 


, , _ 5 ( a o n ) 2a d£ 

P(a 0 S ) - p ■ a 0 4> - p —Tl dt ^ tan 1 ' 

n o a o 


(38) 


(5(a n) 9 d >- 

^(a 0 T) = R P * a 0 * = + — — ~ 2 dt t? tan I . 

n o a o 


(39) 


Designating by a the angular distance of the ascending node of the mean orbital plane from the 
x-axis of the almost ideal system of coordinates we have: 

cr = n 0 at + cr 0 (40) 
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where or is a constant, and taking Equation 27 into consideration 


g - pc os (v-cr) - ( 1 + v) p cos ( f + Cl}) , 

r) - psin(v-o-) - ( 1 + v~) p sin ( f + co~) , 


where we set 


w - K _a o) + n o (y _ a ) 1 • 


Substituting Equations 41 and 42 into Equations 38 and 39 we obtain: 


, N da o n 2 a(l + v) d (^ a o) 

P ’ l a o S ) = P ~J-p FTT 3t Ptanl 


P ' 


( a o T ) 


aa o n 2a(l+ V) d (^ a o) 

W 4 ^ ar-qt-i. 


where we set 

p = -^ COS (T + 03 ) , 

p — 

q - — sin (f + o>) • 

Making use of Equations 44 and 45, the Equations 35 through 37 take the form: 


(41) 

(42) 


(43) 


(44) 


(45) 


(46) 

(47) 


_d 

dt h 


n 0 d ( a o n ) 2(l + ^)aq d U /a o) 

' ‘ “dt tanl , 


/ 1_ e o 


d f 


/” 1_e o 1 2 


(48) 


dT 

ar " n o + 


2 n. 




o 

~P 


1 . _L_\ 

2 1 - e 2 

o A / 


cos f + 


'0 h 2 


1 - e n 


h o 2 


< 9 (a o n) 


5(a 0 n) 


+ ~p sinf • P 5p~ 


4a 




1 + V h 2 ! - ^ 

Sin 03 + — - • —J • q^cos t + e Q j 

1 e o h o 


d(£/a 0 ) 

dt tan 1^-, 
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dW 

dt 


T + 2 - 



+ 


2 n 


o 


y 1 - e o 2 





sin f 


d ( a 0 n ) 

3f 


12 

P 


3 ( a 0 n ) 

cosf-p dp 


+ 



COS Cl) + 


1 + v h 2 

1 ” g 2 h 2 
1 c 0 0 


* q sin 


T ) 


dt 



(50) 


Making use of 


an 

7f 






and of the basic relation 


— sin f 
a o a 




where l is the disturbed mean anomaly 


l 


n„ z + c„ 


we can transform Equations 48 through 50 to a form resembling that contained in the author's 
previous work (Musen, 1963), but with some additional terms. They appear because we choose the 
almost ideal system to be the basic reference frame, instead of the ideal one. We obtain: 


dT 

dt 


d{a o n) <?(a 0 n ) d £/ a o 

+ n 0 y'P + M 1 — + N lP dp + Kj dt 


d»P 

dt 


T+ 2 tt 
0 


1 - e 0 / 


+ 


a(a 0 n) 


d[a 0 n) 


-dZ~~ + n 2 p ~d y 


+ K, 


d£/a 

dt 


d h o d ( a o n ) d ( a 0 n ) d £ /a o 

dt h ” d-l + ^3 P dp + ^3 dt 


where 



(51) 


(52) 

(53) 


(54) 


12 



r 



N 


3 


1 “ e 


2 

0 


e Q sin f 


(61) 


K 


3 


2a ( 1 + v ) q 

A _ e o 2 


tan I 


(62) 


EXPANSION OF THE DISTURBING FUNCTION AND THE FORCE-COMPONENTS 

Let the x-axis of the inertial system be directed along the line of apsides, toward the pericenter 
of the sun, the y-axis be normal to the line of apsides, in the plane of the solar orbit. The z-axis 
will be normal to this plane. Let *>' be the angular distance of the pericenter of the sun from the 
ascending node of the mean orbital plane of the satellite. We have 


- "o' - n 0 yt . 

The position vector p' of the sun can be represented as 


(63) 


r' - p' - P J o'cos(f'+a>') + kxP/o'sin(f'+&j') 


(64) 
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where f' is the true anomaly of the sun. Substituting 


k - Q sin I + R COS I 


into the last equation we deduce 


p' - Pp'cos(f'+to') + Qp' cos I sin(f ' +oi') - Rp'sinIsin(f'+to') . 


The position vector of the satellite can be represented in accordance with Equations 41 and 42 


r P p cos ( f + to) + Qpsin(f+ to) + R £ . 


Making use of the two last equations we deduce: 


r • p' = p p' cos H - p' $ sin ( f ' + to' ) sinl 


and 


A 2 - | r — p' | 2 - (p 2 +p' 2 -2pp' cos H ) + 2p'£sin(f' + to' ) sin I + £ 2 , 


where 


cos H - cos ( f + to) cos ( f ' + to' ) + sin ( f + to) sin ( f '+ to' ) cos I . 


Setting 


_P 

p‘ 


\ 2 


- - S , 


— - cos ( f ' + to' ) 
p 


— sin ( f '+ to' ) - q' 


we can rewrite Equation 67 as 


A 2 = p' 2 A Q 2 


1 + 2 


(^ a o ^ q' sin I 

W a ° s ) 7T 


'5/a 0 


(65) 
as: 

(66) 

(67) 

(68) 

(69) 

(70) 
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From the last equation we obtain: 


~K 



Taking into account 


( 71 ) 


[j/2] 

Pj (x) = J ^ (-l) k A. k xi- 2k . 


where 


A 


i.k 


_j_ ( 2 j - 2k ) ! 

2 j ' k! (j -k)! (j - 2k)! ’ 


we have 


a 


o 


■A 



[j/2] 

('l) k Aj, k (q* sinl)i- 21< 

k=0 


s '2k+l 

A 0 2 i " 2k+1 


( 72 ) 


Making use of the expansion 



L 



(cos H) 


in terms of Gegenbauer polynomials, we deduce 


[j/d 


o s ' " J Z] ( " 1)k . A i- k(q ' sinI)i " 2ks ' 2k ' if)" c n i_k+I/2 ( cosH ) • (73) 

j= 0 n = 0 k = 0 


We set 


P a 0 — 

S. - cos(f-f' +OJ-CJ 1 ') 

1 a o p‘ 


pp' + qq' 


(74) 


15 


p 0 _ 

S 0 - — • — — COS(f + f' + 00 + 0 )') 

2 a 0 p 


pp - qq 


(75) 


S Q - Sj cos 2 ~2 + S 2 sin 2 — 


p “0 


2 a 


o p 


cos H 


(76) 


where the first term in Equation 76 is the main part, 


K 


0 



(77) 


and 




(cos H) 



(78) 


We deduce from Equation 73 


i ~ 0 n- 0 


+ v) n n. 


(79) 


where 


[j/d 


n. 


(-l) k A. k (q' sin I)j' 2k s' 2k T n j ' k + 1/2 


(80) 


the useless term n o 0 = +1 can be omitted. From 


C n X (cosH) 


[n/ 2 ] 


m=0 


(~l) m COS 

N / n . m 


where 


2 n ~ 2m (A, n~ m) 
m ! ( n - 2m ) ! 
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and taking Equations 76 and 77 into account, we deduce the formula which is convenient for obtain- 
ing the expansion of on the electronic computer: 


T 


n 


L 


m = 0 


(-l) m B 


X 

n ,m 


q n - 2m 
^0 


K 


2m 

0 


( 81 ) 


We have the following relations between the series which are a paraphrase of the recurrent 
relations between Gegenbauer polynomials: 


2X + 2n + 2 x 2X+n , x 

n + 2 n+i n + 2 "o *11 


2X 


- (S 0 T^V - K 0 2 T n _ 2 ' ) 


2X 


/ 'T* A.+ 1 _ Q >T* A.+ 1 \ 

n + 2A. \ n ^0 n-1 / * 


n * * ^ „ T \ _ ^ (k 2 - Q 2\t^1 

n S 0 T n-l n V K 0 ) n-2 


(82) 

(83) 

(84) 

(85) 


The application of the relations in Equations 82 through 85 can simplify the actual expanding of T^ 
into trigonometric series with four basic arguments: l, V , co, co' . The expansion on the computer 
can begin with the expansion of the basic expressions 


^ 0 


P T 
cos f 
a o 


Vo 


sin f 


and 




COS f ' 

p 


s in f ‘ 


into trigonometric series in terms of l and respectively of V . Then 

p = coso) - 7} 0 sinco , q = £ 0 sinw + r/ 0 cosa , 

p' = cos a 1 - ? 7 0 ' sinco' , q' = sinco' + cos co' 

are formed. After that in succession are formed 


s 2 . 
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and powers of q' sin I and s' and, finally, fi. n and a 0 / A. The expansions are trigonometric series 
in four basic arguments and with purely numerical coefficients. The numerical values of the mean 
eccentricity of the satellite, of the eccentricity of the solar orbit, of the mutual inclination I and of 
the parallactic factor a 0 can be substituted from the outset. 

For the indirect term of the disturbing function we have 


r ' P 

J ° p ' 3 


p cos H £ 

i ~2 s; * n (f +co )sinl 


or, in a slightly different form, 


a 0 ^3 = a 0 2s ' ( 1 + l/ ) n 0,l- ip n i,0 • 


Substituting Equations 79 and 86 into Equation 16 we obtain 


( 86 ) 


a o n 


1 - 


+ E 


4co . 4 CO 

j = 0 n=0 


v>" n. 

' !■ n 


where 



( 1 + v) 2j 1 a>. 



(1 + v) 2 p 2 



(87) 


(-l) j (2j - 1) ! ! _ 2j _, 

’ ; i s 21 1 

2 i J • 


and for a 0 V we obtain, after an easy transformation: 


a, V = 1 

' m + E 




+ z/) n fi. 

7 J > n 


i l 


2^' s 




( i + uy 2 ^ 1 co. 


a2 D fe) 


2 ^ 

tan 2 I-(l + v)q — tan I 


( 88 ) 
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In the expansions in Equations 87 and 88 we set: 


0,0 


= n 


1,0 


= n, 


0, 1 


0 . 


The expansions to be substituted into Equations 51 through 53 and Equation 18 are 


a(a 0 n) 


al 


j -0 n - 0 


and 


• Efer £ 


2 [a 

a L 2 


<5pj 


dq £ 


f (1 + v) 2 tan" I " (1 + *) Jt • T ' tan 1 


(89) 


3(a 0 n) 


1 '^Ti) m2s 'L ( ‘ 1) 

j=0 


/ \ +” 

(#j 2Z n < i_2 ( itv ) nfi i- n 



n- 1 


(2n + 1) (1 + 12 ) 2n ~ 


OJ 

n 


+ a 


(l + ^) 2 q 2 


tan 2 I~(l+v)q — t an I 


(90) 


d ( a o V ) 

d £/a„ 


1 - 


£V _I V" 


m ' + E 


)m 2 s' (irj ) . a 0 " + i' 2 (l + v) n a n 


i=l 


0 / L-. 

n- 0 


3 v 


L(i + vy (i +v 


j=2 


+ v)‘ 2j_1 OJ. 

1 J 


+ a' 


_ a o 


tan 2 I - ( 1 + v) q tan Ij . (91) 


In fact, only few powers of rJa Q are to be kept in the expansions in Equations 89 through 91. The 
explicit form of fij , n and of T n x which are sufficient for the purpose of the actual expansion up to 
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the fourth power in £/a 0 are 


n i,„ = (q'sinl)T^, 

n 2 ,n = 4(q'sinI)*T n ^ , 

n 3,n = "I (q' sin I) 3 T n 7/2 - (q' sin I) s ' 2 T n 5/2 , 

n 4,n = ^ (q' sin I) 4 T n 9/2 - ^ (q' sinl) 2 s' 2 T n 7/2 + | 


and 


•T 1/2 
X J 


S 0 ■ 


T 1/2 
2 


1 5 2 „ 1 v2 

2 2 K o 


T 1/2 

1 3 


5,3 

T S o 3 - 2 s oV 


35 . 15 „ , 3 . 

-8 V - 4~ S o K o 2 + 8 K o 


T 3/2 

1 


3S 0 , 


T 3/2 
2 


c 2 . 1 p 2 
2 b 0 2 K 0 


3/2 


35 


15 


T S o 3 * T S o K 0 2 


5/2 = 


5S n 


35 , 5 , 

T S o 2 ' 2 K o 2 


7/2 = 


7S 0 , 


s' 4 T 5/2 

n 
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THE PROBLEM OF INTEGRATION 


The integration of the differential Equations 51 through 53 and Equation 18 can be accomplished 
by their reduction to integral equations and then by applying the method of iteration. The operator 


where 


D 


_d 

dt 


d n o z a , d . d d 
dt + n ° aV + “ a* + “ ao,' 


“ - n 0 (y _ a ) 


co' - n # asecl , 


can be represented as a sum of the operators 

D □ D 0 + D t 


(92) 


where 


and 




. d 

+ CO ■+ CO 

d co 


d 

dco ' 


(93) 


dn o Sz <9 

dt el (94) 


and the value of d n Q 8z/dt is given by Equation 97. Each of the differential Equations 51 through 53 
and the equation for d n 0 Sz/dt have the form 

-t* , CO, co' ; 1 + V, 5 / a„; h/h 0 , T, W) . 

They can be re-written as integral equations of the form 


<P - D 0 1 (<D-Dj <p) 

to which the method of iteration can be applied. 

The process of formal integration of Equations 51 through 53 and of the equation for the per- 
turbations of the mean anomaly is thus reduced to the simple application of the inverse of a linear 
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partial differential operator with constant coefficients to trigonometrical series. In the lunar 
problem each series is either a purely sine or a purely cosine series in four basic arguments. If 
the result of the integration is a cosine series, then an additive constant of integration must be in- 
troduced. Such constants appear in h 0 /h and in T. The series for w andn 0 Sz are pure sine series 
and additive constants are equal to zero. The value of n 0 y at each integration step is obtained in 
such a way that no constant term appears in the right side of Equation 52. This is done to avoid 
the presence of a secular term in <P. The iterative process is continued until the computation of 
series of 1 + v, n 0 Sz , £/a 0 , h 0 /h, T, f, y and a leads, with accepted numerical accuracy, to the same 
results. 

At this step we shall introduce Hansen’s w function in a manner similar to the classical theory. 
We set 


w 



cos f + 


5 o) 


p — 

W — sin f . 
a o 


(95) 


The analogy is however, only external. The w function as defined above is not identical with the 
corresponding function of the classical theory. The classical W is associated with the perturbations 
of the auxiliary satellite moving in the osculating orbital plane and the influence of the periodic 
perturbations appears in W only indirectly. Our w is associated with the perturbations of the auxil- 
iary satellite moving in the mean orbital plane and the influence of the periodic perturbations of 
the osculating orbital plane is reflected in our W directly by the presence of terms of the form 
K i (d £/a 0 /dt ) in Equations 51 through 53. 

Thus, in the present theory the W function represents a more or less formal device whose in- 
troduction is suggested by the form in Equation 8 of the equation of motion of the satellite in an 
almost ideal system of coordinates. Equations of type 8 permit one to transfer to our problem 
all the formalism of the classical theory, or at least that which concerns the perturbations of the 
radius vector and of the mean anomaly of the projection of the satellite on its mean orbital plane. 
We have, for example, the Hill formula (1881) in the form 


or 


Let 


d n Q Sz 
dt 


W + I' 2 

1 


n o y P 2 
a o 2 


d n Q 82 
dt 


n 0 W + 


n Q 7 (1 + W) 

71 .,2 





[T] , [<P] 


(96) 


( 97 ) 
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be the values obtained by the formal integration toward the end of an iteration step. Of course, 
these values will be different for each specific cycle of the iterative process, but in order not to 
complicate the writing we shall use these symbols indisc riminantly for all cycles. Similarly as in 
the previous modification of Hansen's theory and in accordance with the remarks made above, we 
have 


"0 

“FT = 1 + c 


’V 
1 + L h_ 


= 1 + A 


T = c 2 + [T] and <p = [<p] 


where c and c 2 are the constants of integration. From Equation 30 we obtain 

H = (- 3Cj - § e Q c 2 ) + [~] , 

where we set 


(98) 

(99) 


(100) 


rh 


[~] 


? 0 [T] + 2(A 2 - A 3 + • • •) 


( 101 ) 


and A 2 , A 3 , • • • may be taken from the previous iteration. Taking into account Equations 99 and 100 
and setting 


[w] 


cos f + I e 0 ) 


+ |^- cos ft o e„) [T] + sin f M 

d o 


we have 


/, , 3 ^ 

( P - 3 \ 

\3 c i + 2 e o c 2 1 

1 + n o c 2 cos cp + 2 e 0 j 


+ [w] 


and Equation 97 takes the form 


d n Q §2 
dt 


- n o ( 3c i + | e 0 c 2 ) + n o c 2 cos cp+ § e Q ) + Q 


(102) 


(103) 


(104) 


where 


Q 


[w] + 


n 0 v2 C 1 +w ) n 0 y 


A 


- e 2 
c o 


(105) 
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II 


The constants Cj and c 2 must be determined in such a way that the constant term and the term of 
the form A cos l do not appear in the right side of Equation 104. The determination of these con- 
stants is fully discussed in the author's previous paper and therefore is omitted here. The integral 
equation to be solved by the method of iteration has the form: 


or 


n 0 Sz 


f„( 


3c, + 


l) +n o c 2 (‘r 0 cos cp+ 4 e 0 ) 


+ Q~ D i n 0 Sz 


} 


(106) 


n o 52 


/ / 3 ' 

1 ( P - 3 > 

dn Q Sz (9n 0 Szl 

\~ n o l 3c i + 2 e o C 2 > 

1 +n o c 2{-^ cos q>+ 2 e o ; 

^ +Q dt dl J 


(106') 


The values of v 2 , w and of n 0 Sz in the right sides of Equations 105 and 106 again can be taken from 
the previous cycle of iteration. 

The perturbations in the radius vector are obtained without integration, as in the author's pre- 
vious work, using the formula 


2 ( A ' 


W) 


(A +W)y 


The integration of the differential Equation 18 requires some special considerations. We re- 
write this equation in the form which favors the application of the variation of constants: 


D o 2 


u. 

/ a o 


= Z 


(107) 


where 


and 


Z 


„ 3 a o V d(l+v)p 

n o d 1/ a Q + 2n o a dt tan 


- 



D 2 = D 0 D 1 + D 1 D 0 + D l 2 • 

The solutions of the homogeneous equation corresponding to Equation 107 are 

£i l 2 


(108) 
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By applying the method of variation of constants we obtain the integral equation 


i ( 109 ) 

which we shall solve, by the method of iteration. At each iteration step the approximation to n 0 a 
must be determined in such a way that no term of the form A sin (f. + &>) appears in 4/a 0 . The con- 
vergence of the iteration process for xj a 0 , because of its smallness, is fast as compared to the 
speed of convergence for n 0 Sz. In fact, only a few powers of £/a 0 are to be kept in the right side 
of Equation 109. 

COMPUTATION OF THE POSITION VECTOR 

The computation of the position vector at the given time t is based on the numerical evaluation 
of 1 + v, n 0 8z and c using the corresponding series. It starts with the determination of l and n 0 Sz 
by means of successive approximation using the expression 

1 -= n 0 t + c 0 + n 0 Sz ( 110 ) 


and the series for n 0 Sz, and beginning with 

1 = n 0 t + c 0 


as the first approximation. The numerical values of three other basic arguments, l, m and co' are 
obtained from Equation 92. After l and n 0 Sz are computed we evaluate the numerical values of 
1 + v and 4 from their series. Then the position vector of the satellite is obtained by using the 
formulas: 


A 3 (-&)') ‘ Aj (cos I) • A 3 (o>) 


(1 +l O a o ( cos E ~ e o) 
(1 + v ) a 0 y'l - e 0 2 sin E 

4 


E-e Q sinE - -£ , 


where A t and A 3 are matrices of the form: 


+1 

0 

0 



+ cos a 

- sin a 

0 

0 

+ cos a 

- sin a 

A 3 (a) 

= 

+ sin a 

+ cos a 

0 

0 

+ sin a 

- cos a 



0 

0 

+1 


l 

, 
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There is no real necessity to invert the series for 1 + v , n 0 S z and £ in order to represent the per- 
turbations in terms of four arguments linear in time, because the solution of Equation 110 for l 
can be accomplished very rapidly and without any inversion using an electronic computer. 


CONCLUSION 

By referring the motion of the satellite to an almost ideal system of coordinates it is possible 
to put Hansen lunar theory into a form which gives the coordinates in a more direct manner. The 
number of auxiliary constants is less than in the classical version and it is also less than in the 
author's previous version. The necessity of a literal expansion of the disturbing function is re- 
moved and the development is made a purely numerical one. 

We decided to use the general perturbations method to represent the motion of satellites be- 
cause such a method, when it can be used, gives a better physical description of the motion as an 
oscillatory process, as compared to the method of numerical integration. The proximity of Jupi- 
ter X in 1967 to its predicted position serves as a check of the theoretical thought and shows that 
the machines can be used to perform analytical expansions in celestial mechanics. 


Goddard Space Flight Center 
National Aeronautics and Space Administration 
Greenbelt, Maryland, December, 1967 
188 - 43 - 01 - 01-51 
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